一次元圧密の熱移動解析のコードです！

chatgptを参考にしました。

import numpy as np

import matplotlib.pyplot as plt

# 土の物理パラメータ

k = 1.0e-9  # 土の透水係数 [m/s]

mv = 0.8e-4  # 体積圧縮率 [1/Pa]

es = 1.0e-3  # 土の初期間隙率

# 熱パラメータ

rho = 2000  # 土の密度 [kg/m3]

c = 840  # 土の比熱 [J/kgK]

kappa = 1.0e-6  # 熱拡散率 [m^2/s]

# 初期条件

initial_temperature = 15.0  # 初期温度 [C]

# シミュレーション設定

L = 10.0  # 解析領域の長さ [m]

dx = 0.1  # 空間ステップ [m]

dt = 10.0  # 時間ステップ [s]

N = int(L/dx)  # 空間分割数

T = 86400  # シミュレーション時間 [s]

steps = int(T/dt)  # 時間ステップ数

# 空間と時間の配列初期化

p = np.zeros(N)  # 圧力配列

T = np.zeros(N) + initial_temperature  # 温度配列

def update_pressure_and_temperature(p, T):

    new_p = np.copy(p)

    new_T = np.copy(T)

    for i in range(1, N-1):

        # 圧密方程式の差分形式

        new_p[i] = p[i] + dt * (k/mv) * (p[i+1] – 2*p[i] + p[i-1]) / (dx**2)

        # 熱伝導方程式の差分形式

        new_T[i] = T[i] + dt * kappa * (T[i+1] – 2*T[i] + T[i-1]) / (dx**2)

    return new_p, new_T

for step in range(steps):

    p, T = update_pressure_and_temperature(p, T)

    if step % 100 == 0:

        plt.figure(figsize=(10, 5))

        plt.subplot(121)

        plt.plot(np.linspace(0, L, N), p)

        plt.title(f”Pressure at time {step*dt} s”)

        plt.xlabel(“Depth [m]”)

        plt.ylabel(“Pressure [Pa]”)

        plt.subplot(122)

        plt.plot(np.linspace(0, L, N), T)

        plt.title(f”Temperature at time {step*dt} s”)

        plt.xlabel(“Depth [m]”)

        plt.ylabel(“Temperature [C]”)

        plt.show()

例えば，以下のような📉が出来上がると思います。

パラメータを変更してシミュレーションしてみてください！！！